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I The post-Newtonian (pN) description of the dynamics of compact binary systems is the topic of actual research. 
(2jr)' Due to the strong nonspherically symmetric gravitational interactions, those objects are supposed to be "secure" 
^ ' sources for the detection of gravitational waves. Currently, LIGO, VIRGO, and GEO600 search for the last seconds 
or minutes in the life of those sources when the gravitational wave (GW) emission frequency enters the bandwith of 
the mentioned detectors. 

The computer resources - in contrast - are currently unable to create numerical GW templates for the early stage 
of the binary inspiral (this case is important for LISA) , where hundreds or even thousands of GW cycles have to be 
simulated which makes it necessary to propose an analytical prescription of the orbital evolution, and in the case of 
^ spinning compact binaries, the spin evolution as well. An equally essential ingredient of the GW data analysis is the 
I ' transformation of the GW signal into the frequency domain. There exist numerous more or less optimized numerical 
^ . routines to convert the time domain signal into the frequency or Fourier domain. To economize on computer resources, 
it is also reasonable and desirable to provide analytical Fourier domain wave forms for the data analysis community. 
Numerous authors have investigated the performance of circular inspiral templates and their analytical Fourier-domain 
' pendant, which have been set up to a certain standard in the literature thanks to the work of Daniour, Blanchet and 
I Iyer to compute higher pN order corrections to the GW energy loss [H-Ql • 

i This work is inspired by Q and @ and is a direct sequel of Q. Galt'sov et al. provided a formal frequency-domain 
' decomposition in multiples of two irreducible frequencies: the radial frequency fr, describing the elapsed time form 
. one periastron passage to the next, and the frequency associated to the periastron advance parameter, f^, but used 
this decomposition only for the computation of the far-zone energy flux, where some mathematical aspects made it 
unnecessary to compute the Fourier coefficients directly. 

Ref. provided numerical insights to the case of no radiation reaction (RR) with the use of essentially the same 
frequencies that appeared in 

Yunes et al. @ established analytic Fourier-domain inspiral templates for eccentric binary emission for the exem- 
plary case of Newtonian equations of motion (EDM) and the leading-order GW amplitude. Those authors used the 
, expansion of the Kepler equation and related trigonometric function combinations of the eccentric anomaly in terms of 
^ ' Bessel functions and applied this to a series expansion of the leading-order GW amplitude in the orbital eccentricity. 

What we like to do in this article is the generalisation of Yunes et al. to IpN order of the conservative orbital 
dynamics. The work will include, in addition to Yunes et al., the effect of periastron advance and also the GW 
amplitude corrections to IpN order relative to the leading order quadrupole approximation. We present the results 
in terms of irreducible positive frequencies fr and f^. Afterwards, we compute the Fourier domain for the case when 
RR is taken into account. 

The paper is organized in the following way. Section U gives an overview over the solution to the IpN-accurate 
orbital motion. Section |TT] summarizes the transverse-traceless (TT) projection of the far-zone GW field up to IpN 
order corrections to the leading-order quadrupole field. The vital components of the Fourier decomposition of the 
GW field are summarised in Section IIIIl Section IIVI depicts how to decompose the terms of the multipole moments 
of equal structure. The Fourier domain multipoles themselves, incorporating purely IpN conservative dynamics, are 
given in |Vl The effect of RR is incorporated in Section I VII using the method of the stationary phase approximation 
(SPA). 
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I. ORBITAL MOTION OF COMPACT BINARIES 



Before we start computing the GW forms, we give an overview on the sohition to the orbital EOM in a quasi- 
Keplerian parameterization (QKP). The QKP for non-spinning compact binaries up to and including 3pN point 
particle (PP) contributions can be found in @ and references therein. We restrict ourselves to the IpN part thereof. 
The QKP for nonspinning compact binaries to IpN is the following: 
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In the above Eqs. ([ij-lll]), e counts the inverse order of the speed of light c, = c~^, r is the scaled relative separation 
of the two bound objects, cj) is the orbital phase, v is the true anomaly, k is the periastron advance parameter, and et 
are some "radial" and "time" eccentricity. The latter appears in the Kepler equation (KE) , Eq. ^ , connecting the mean 
anomaly M (which is directly proportional to the elapsed time t with the proportionality factor n = 27: / = 27r/r, 
where is the time between two consecutive periastron passages) to the eccentric anomaly m, see Q for geometrical 
insight. The shorthands 77, and \E\ denote the symmetric mass ration rj = {17111712) / {mi + 7712)^ ("^i a-nd m2 are 
the individual masses with = mi + 7712); the reduced mass, and the absolute value of the orbital binding energy. 

This parameterization is essential for the computation of the higher-order time derivatives of the radiative multipole 
moments. We can easily obtain them using the KE. How the orbital parameterization effects the far-zone gravitational 
field will be depicted in the next sections. 



II. THE IPN-ACCURATE TRANSVERSE AND TRACELESS FAR-ZONE RADIATION FIELD FOR 

ECCENTRIC COMPACT BINARIES 

We start collecting the expressions of the transverse-traceless far-zone field h'^"^ in Junker and Schafer They are 
given in tensor-spherical harmonics and it is remarkable to note their structure as they appear as ij^^ = e"""*^ .fim{u{t)) 
and sill = e"*""^ 9im{u{t)). Taking as basis 
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where R is the distance from the observer to the source, T^^^''™ and T^^^''™ are the "electric" and "magnetic" tensor 
spherical harmonics [l^, the superscript {i) allowing i = 2,3,4 denotes the i*^ time differentiation. They read - up 
to the required overall order of = cr^ in Eq. PU]) - 
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where is the reduced mass and A{u) := 1 — 6(C0Su. The above multipole moments contain terms with certain 
equal structure. Evaluating the exponential functions e"™**^ with m = 1,2,3,4 to IpN accuracy, one obtains factors 
of the form A{u)^'^ and sin(u)/^(u)". ^ Those terms can be expanded first in sin and cos series of u and afterwards 
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in terms of M with the help of Eqs. (jTB]) and (|18p . We wih calculate the expansion coefficients in the next section. 
Some calculation steps will be required repeatedly and we will collect them as ingredients for a frequency-domain GW 
recipe. 

III. INGREDIENTS FOR FOURIER DECOMPOSITION 

The Fourier-domain GW form requires computation of relatively simple structures which will combine in the tensor- 
spherical harmonics. Let us start with the most fundamental quantities which we compute from the scratch on the 
one hand and some which we collect from the literature on the other. 

1. The inverse scaled radial separation with an arbitrary integer exponent n > 0: 
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Equation (jl2p is proven in Appendix |X1 Note that from here onwards, e can be set or et, depending on the 
context. 

2. Trigonometric functions of multiples of the eccentric anomaly u in terms of Bessel functions (see [llj, p. 555), 
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Note that the above Eqs. (|T6|) - ([20)) are only valid taking et because of the structure of the KE. 
3. The decomposition of v in M (see [l^l, P- 33), 
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4. The decomposition of the IpN accurate Kepler equation itself in terms of the mean anomaly (see p. 553), 
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which converges for all < 1. These inputs are sufficient to obtain a IpN-accurate F-domain decomposition of 
the IpN far-zone GW field including f pN-accurate orbital dynamics. More complicated expressions will be based on 
those above and can be computed more or less laboriously. A detailed appendix of this article will give some proves 
for resummation rules for these terms. The next section will deal with certain exponential functions including the 
IpN-accurate phase. They will be combined with the fast oscillatory terms in and g;„i and have individual Fourier 
decompositions as shown below. 

IV. DECOMPOSITION OF THE ELEMENTS APPEARING IN THE MULTIPOLES 



This section is subject to the Fourier decomposition of single terms in the multipole field expressions. We will apply 
the aforementioned steps to obtain the individual Fourier decompositions. The "magnetic" number (the factor of —i (/) 
in the phase exponentials of Eqs. (|llap - (|llsp ) which we shall call m in the multipoles dictates terms of the form 
g-im0^ which have to be expanded in e. 



2 arctan 



:= 2 arctan 



1 + et u 
tan — 

1 - et 2 




-im(pa „-«"i(l+k)fe^ 



-irrKpQ 



-im(po 



1 — 60 cos U 

COSU — 60 
1 — 60 cos U 

COS u — 

1 — 60 COSU 



e 

sinw 



1 — e^cosu 



~ smu 

1 — e~i 

''^ 1 — e^cosu 



svau 



^4> 1 



60 < 



(26) 
(27) 



exp 



ik\M + ^ Gj{et) sinjM 



exp [— imkM] 1 — imk Gj{et) sin jM 



(28) 



For the first line we took the definition of w as a function of u and et, Eq. (|4]), and Eq. ([3]). In the last line above we 
truncated the exponential of the sin series after the linear order in k, because k is already of IpN order and used the 
Fourier-domain expansion Eq. ([55]) for v. We are allowed to do this because there is no secular contribution coming 
up in the exponent due to its periodicity of 27r in M and its average to zero in the interval [0, 27r]. 

In the multipole moments the numbers m = (0, 1,2,3,4) appear, so we preserve all contributions with these m up 
to the order of e required. 
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Those terms will mix with /;„ and gi^. Let us get more precise now and start decomposing the A{u) " term. With 
the help of Eq. we get 
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where we have collected those terms with the same positive frequency. Appendix |A] will give deeper explanation how 
to select the frequencies with the same absolute value. The same contribution - being multiplied with sin u - will be 
decomposed analogously, but having sinmw this time as an odd function of u. 
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Remembering Eqs. ([16]) and (|T8l) . we get 
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The abbreviations and stand for the j"^ contribution of A(w) "and A{u) " sinu. Of course, both expressions 
wiU not take their place alone, but will be multiplied with the expansion of the e*"*"^"* term due to Eq. ([25]) . We 
Fourier decompose the following products of series as they appear in the field multipole moments. The first quantity 
consists of two pure sin series. Factor 1 is the expansion of sin(w) A(w)~" in M and factor 2 is the non-secular part 
of the exponential of the w(/)-term. 
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The second quantity consists of a sin and a cos series. Factor 1 is the decomposition of A{u) " in M and factor 2 
is - again - the non-secular part of the exponential of the w(Z)-term. 
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The two above summation formulas will also be proven in Appendix |D] and [C] Please keep in mind that the (n) 
is always present in those terms and reminds the user of keeping the right exponential of and that those 

functions always depend on the value of e. Additionally note that the statements "for j > 1'' and "for k > f ensure 
that, for example, the first term in (4.19) is not able to contribute to the first harmonic in / and that the indices of 
the decomposition coefficients do not become zero, or else, these elements are set to zero automatically. 

The next step is to sum up and summarise the terms appearing in the multipole moments. Using Eq. ()22[) . they 
also contain absolute static elements and will be listed below. 
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We are now in the lucky position to decompose all point particle contributions to the radiation field amplitude in 
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and use them from now on in the multipoles, where n and q will also be allowed to be 0. We finally rewrite our wave 
forms in the following manner, 
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with some K^rn] and Krimi to be determined. 



V. DECOMPOSITION OF THE MULTIPOLE MOMENTS 



We will decompose the multipoles in terms of the functions J^[nq] (u) introduced above next. These functions will 
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/3[3i] := -~et, (80) 



A31] - -^yi-ef, (81) 

"(2) 

520 = 0, (82) 

"(3} / o_ \ 4 f ^ ^ ^ 

/ 33 = 8 - 77(mi - m2)(-i?)3/2e-^^''^e-'3'^" ^ 7[.3]-^[fe3](^i) + E 7[.3]-F5[fc3] («) , (83) 



,A:=0 fc=l 



7[03] := '-^ (84) 



7[i3] := ^ (85) 



^(l2-7e2)yT^ ^ ^ 

7[23] ^ ^1^^ (86) 



. ^v/^^(e?-l) (e^ + 4) 
7[33] ■= ^^3 , (87) 



10zVT^(e?-l) 



2 



61^1^ (e?~ 1)3 ^ ^ 

7[53] := '-^ (89) 

7[i3] := 4-3, (90) 

7[23] := 0, (91) 
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7[33] := -^-^ + Y' (92) 

7[43] := ^ * 2 , (93) 

7[53] := - ^*2 ^ ^ (94) 

13) 

/ 32 = , (95) 

(3) FItt \ ^ 1 

/31 = e-'''^e-^^°8^/ — (~i?)3/2(mi-m2)r;x <jE7[fei]J-[,i](^) + E7[fei]J-s[fei](7.)|>, (96) 



. /c=0 fc=l 



7[o, - (97) 



7[ii] •= : (98) 



11 



7[21] 


^ ■ /-I 9 

bet 


(99) 


7[31] 




(100) 


7[ii] 


:= -1, 


(101) 


7l211 


:= 0, 


(102) 


7[31] 


■■= -U^t-i). 

g V I / 


(103) 


(3) 
/ 30 


= 0, 


(104) 


(3) 
533 


= 0, 


(105) 



?32 = e-^■''^^e-^'^''y^E'-,/T^,^,{^-3v)x\J2^[k^^^^^ (106) 



^fe=2 k=3 



'5[22] := (107) 

6* 



^[32] := (2^-|)\/l-e,^ (108) 



<5[42] := ( ^-10?:) Jl-e^ (109) 



6zyr^(e?-l)2 ^ ^ 

'5[52] (110) 



e 



2 

Si32] ■■= et, (111) 

et 

^[42] := 7-4et, (112) 

£•[52] ■= , (llo) 

et 

13) 

531 = 0, (114) 

530 = -16^^^'etVl - etV(l - 377)^[3o]^s[03] (") , (115) 

^[30] := 1, (116) 



6 (et - 86? 



sV(i - 3?;) X I.Y.CmJ'miu) + ^C[fe4i-^s[fe4](u)| , (117) 

lfc=o fc=i J 



C[04] := 4 , (118) 

- "W'y^^) . ,119, 

^t 

C[24l := , (120) 

3e? (ef + 8) - 8 , , 

C[34] := * ^ * 4 19 > (121) 
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280 (e? - 1) 4 

Cm := (124) 

120 (e? - 1) 5 

C[74] := 4^^' (125) 



C[i4] -3 , (126) 

C[24] := 0, (127) 
C[34] := 2zv/T^(9e^^-46e^ + 44) ^ ^^^^^ 



8iJl - (5e4 - 17e^ + 12) , ^ 

C[44] := ^^-1 ^ ^, (129) 

12z (1 - e2) 5/2 (e2 - 8) , , 

C[54] ^ '-y^ (130) 

~ 160i (1 - e?) 7/2 

C[64] 3 , (131) 



e 



- 120z (1 - e?) ^ ^ 

C[74] := 3 , (132) 



~^ ^ 

/ 43 = , (133) 

(4) 8 f 1 

/42 = e-'2kMg-i2^o^V2^£;2^(l_3^)x ^C[fc2]-F[;c2](^)+EC[fe2]-Fs[fe2](^i) , (134) 

lfc=0 fc=l J 

12 

C[02l 6 2' (135) 

C[i2] := ^-6, (136) 
16 

C[22l := --11, (137) 



e 



+ (139, 

12(e,~l)3(e, + l)3 

C[52] = 2 ' (140) 



C[l2] •- , (141) 

et 

C[22] := 0, (142) 



» (3e, - 4) (3et + 4) 

C[32] •= , (143) 

et 

8^(l-e?)3/2 

C[42] := ^ '-^ , 144 

et 

8i (1 - e?) 3/2 

C[52] := ^ 145 

et 

/ 41 = , (146) 

?40 = ^yf^V(l - 3r/) |5^C[fc2]-^[fc2](") J ' (147) 
C[oo] := 6, (148) 
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C[io] 

C[20] 
C[30] 



-6, 
-5, 

5 - Se? , 



(149) 
(150) 
(151) 



It seems that there are divergences for et — ^ 0, but the sum of the divergent terms is finite. Each decomposition 
term will contribute to sinjM and to cosjM, allowing for j = [0,oo), with its own values [nm\ in boxy brackets 
for the multipole moment having the same special "magnetic" number m and n labeling contributions with A(u)^". 
Using this and, not less importantly, Eqs. and ([50]) to separate cxp(— imkM) from the products of A{u)^" and 
sinu A{u)^" with cxp(— jmwet ), the moments have the Fourier decomposition 



(2) 
/ 22 



g-2i0Og-2ikM 



27r 



(152) 



\[fc,2] 



Kj=l 



( 5 



fc=l 
5 



Y COS jM <j Y «[fc2] {Mi + "[fe2] iSc) 
J=0 



[fc,2] 



, A;=0 



fc=l 



(2) 
I 21 



0, 



(153) 



(2) 
/ 20 



^' 15 ' 



^sinjMi ^a[fe2](A 



N[fe,2] 



^cosjMi^a[fc2](A) 



[fe,2] 



\3 = ^ 



.k=0 



j=0 



(154) 



(2) 
S 22 



(155) 



(2) 
5* 21 



g-#Og-ikM 



|y|V^rr;?rK-i?)'/'(mi-m2)x 

(oo r 3 

5]sin.7MK]/3[,i](A)f''' + A3i](5.) 
j=l U=2 

oo r 3 

+ ^cosjMi ^/3[fci](A) 

i=0 Ia;=2 



[1,1] 
j 



/3[31](5c)i- 



(156) 



(2) 
S* 20 



0, 



(157) 



(3) 
/ 33 



g — 3i0og~3i k M 



sJ — i-Ef/^ (mi - ma) r/x 



J^sinjM 5:7.3(A)['=-^1 +^7.3(5.)[^^1 

oo 5 5 

-^cosjM ^7.3(A)['=^^1 +E7'.-3('5c)['='l 

j=0 Ia:=0 fe=l 



(158) 



(3) 
/ 32 



0, 



(159) 



(3) 
/31 



(3) 
/ 30 



g-i</.0g-4kM 



27r \3/2 , 

— (-E) (mi - ma) r/x 
35 



3 



^^sinjM ^7.i(^.)['='^l +E7/ci('5.)['^^l 

[k^O k=l 
Qo f 3 3 

-J^cosjM E7.i(-4c)['^^^l +E7fci(5c)['^^l 
i=o lfc=o fe=i 



0, 



(160) 



(161) 
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(3) 

^33 = 0^ (162) 

"^3^ ZI~'ZZZ, 8 l2n 



S32 = e-2#og-2.kM _^=_E2^l-e^tfi{l-3f])x (163) 

(oo C 5 5 

3 = 1 Vk=2 fe=3 J 

oo r 5 5 ~) \ 

+£cosjmK]4i(A)['='^1 +^4i(4)['=^l , 

j=0 lfc=2 A:=3 J / 

g3i ^ 0^ (164) 

530 = -16 J -^E'^et^l - 6*2^(1 - 377) X 

(^^^sinjM|(j[3o](5,)f +£cosjm|(5[3^^ , (165) 

(4) 4 /27r 

/44 = e-4#0g-4«kM ^ /f![^2^(l_3,^)x (166) 

9 V I 

sinjM/^ a4(A)['^^^l + E Cm(5.)['=^1 1 

Ia:=0 k=l ) 

+f;cosjM|^CM(A)['^-'^l +ECm(5.)[^-4) , 

3=0 \k=0 k=l J / 

(4) 

/43 ^ 0^ (167) 

?42 = e^2»0oe-2zkM Av2^£;2^(i_3^)x (168) 

63 

(ex. C 5 5 

sin y: CkMsf^'^ + E c^^i^sr^ 
3 = 1 [k=0 k=l } 

oo C 5 5 \ 

j=0 lfe=0 fc=l J/ 

(4) 

/40 = 0^ (169) 

^« = ^yf|£;|V(l-3r;) (f] sin jm|e a.o(A)ri | + £ cos jM j ^ Cfco(^c)r | ) • 070) 

\i=l U=0 J j = l U-=o J/ 

The Fourier transformation of the above gravitational waveforms is easily done. We take the Fourier transformation 
of the sinjM and cos jM terms, 



^ J e"""'''^sinjMe*"*dt = i^S{jn + kmn - cu) - i^S{jn - kmn + cu) , (171) 
J e-^^'^eosjMe^^Mt = y^(5(jn + k?7in - w) + y^(5(jn - km7i + w) , (172) 

and apply this to the sums over j. How many terms of the infinite sum combinations are necessary to determine the 
multipoles up to a certain arbitrary order in the eccentricity of the binary system is discussed in detail in Appendix [Fl 
Note that in Eqs. (|153p to p70p . the expansion coefficients are not longer dependent on M. They are functions of et, 
n and the "magnetic" number m and include multiple sums as we keep in mind the previous sections. In IpN-accurate 
orbital dynamics, they are constants. Due to RR, they will become slowly varying functions of time. All the terms, 
where this condition is satisfied and which will appear together with an explicit phase factor, can be treated with the 
help of the SPA. All other ones which do not have a phase factor have to be treated individually. 



15 



VI. THE EFFECT OF RADIATION REACTION AND THE STATIONARY PHASE APPROXIMATION 

OF THE GW FIELD 

A. Radiative dynamics 

Having the conservative evolution under control, we can apply the Fourier decomposition to the radiative evolution. 
This is done by separating the full orbital evolution into two time scales: the orbital and the reactive timescale. The 
latter is - for our calculation - assumed to be much larger than the orbital scale. This is equivalent to the assumption 
that the rate of change of the orbital frequency is small measured over one orbit. The time dependence of the orbital 
frequency is closely related to the loss of energy and orbital angular momentum due to the GW emission. Peters 
and Mathews (Tsl . [3| proposed a relatively simple model of the binary inspiral in the approximation of slow motion 
and weak gravitational interaction for arbitrary eccentricities < e < 1, which has been standard for some time. 
Blanchet, Damour, Iyer, and Thorne made endeavors to obtain higher-order corrections to the quadrupole formula 
P, [13, [la, [l^. The far-zone fluxes of energy and angular momentum to 2pN order can be found in 0. These 
corrections have shown up to be unrenounceable for the data analysis community. For the time being, we will restrict 
ourselves to the Peters-Mathews model for an exemplary calculation. Higher-order corrections can be included in a 
forthcoming publication. 

In the conservative case, n and were constants of motion, but when the orbit shrinks due to RR, both elements 
will follow coupled EOM connected to the loss of orbital energy and angular momentum [13, [3| , 

n = AA(n,et), (173) 
et = £{n,et), (174) 

symbolically. Note that above Eqs. hold after averaging over one orbit. Thus, n will not simply be M = n{t — to), 
but will satisfy 



M{t) = f n{t')dt' , (175) 

J to 

M{t) = n{t). (176) 



The Fourier-domain wave forms will not longer follow Eqs. (|153p - (|170p and will require some more considerations. 
We will apply the SPA at the first order to approximate the frequency-domain wave form first to each single frequency 
separately and will later sum up all the terms of the discrete decomposition as a virtue of the linearity of the Fourier 
integral with respect to the integrand. 



Suppose an integral of the form 



B. SPA of the GW signal 



r+oc 



Then, having found a stationary point t*, where the phase defined as $ := 2 7r/ — (t>{t)t has zero ascent, $(to) — 0, 
with the assumption of the right behavior at the boundary, the integral takes the form 

h{f) « ^/2^4^e-'(*(**)-2-/**+-/4) (178) 



'<P{t*) 

In the case of our GW, the integral turns out to be 



/.(/)_ = e-""^*« / + cosjM) e—"-^ e^^^^* At , (179) 



(") 2 1"^^ 



Kf) 



C,nm 



2tt 



/ + 00 ^ 
-oo 2 



16 



-mi(f>Q 



2tt 



+ 00 



i(jM+27r/t-mkM) 



■i(-jM+27r/t- 



(180) 



where the subscript "C" shall denote the contribution of cosjM. The reader should carefully note that C and S are 
functions of the elapsed time, as the orbit decays due to RR. We assumed that the binary system evolves far away 
from the last stable orbit, such that the periastron advance parameter is much smaller than unity, k ^ 1. Thus, 
nk ^ j with j = 1, 2, 3, .... The case j = will be discussed in one of the upcoming subsections. In the first part of 
Eq. pSOp . therefore, there exist no points of stationary phase for j > and this term vanishes. The second term in 
(|180p . having the exponential argument i (— jM 



2Tr ft — mkM) will contribute, since the phase $,„j defined as 
= 27r/t - jM - mkM 



has a stationary point at t* where 



27r/ — jn — mkn ~ . 



(181) 



(182) 



Here we note that k = ©(e^) and will be consistently neglected. We want to compute t* and its contribution in 
the following lines. The question arises how to find out t* without solving Eqs. (|173p and p74p numerically. It is 
answered in a simple manner. We search for the solution n — n{et, eto)-, ^to being the eccentricity ed t ~ Iq, of 



dn 
det 



n 



insert it into Eq. (|174p . 



et 
et 



£{n{et),et,eto) , 
et{eto,t- to) , 



(183) 



(184) 
(185) 



and invert Eq. (|185p to find (t — to) as a function of et and eto, where for the stationary point the value corresponds 
to the solution of the equation 



= 2TTf - jn{el,eto) - 77ik(?i(et , eto), e^* ) n(et , eto) ■ 



(186) 



We see that and thus t* depends only on the magnetic number to, the summation index j, the eccentricity and 
mean motion at to-, eto and no, and of course the frequency /. The sin term of all those calculations will have the 
same stationary point t*, such that the full Fourier-domain wave form will have the following appearance: 



(n) 

h (/)n,m>0 


' j=0 


\- 5oj)C*^j^ 


1 


e<*".j(C,)-V4) 


— im4)Q 


(187) 




^ 1 


(n) 

h {f)n,m=0 






1 j('foAtL)-^/4) 






(188) 


-V 


/4>Oj{t*o^) 








4'mj {tmj ) 


:= n{j + TOk) \ t=t'^. 












(189) 




Smj\t=tl^j , 












(190) 


/-^* 


■ = C,nj\t=f^^. ■ 












(191) 



The addend 5oj appears only at one special place in Eq. (|187p because in Eq. ()180p both terms, the first and the 
second, contribute to the cos terms for j ~Q and the sin terms cancel each other. For computing the energy E which 
is a usual prefactor with some exponent in all the multipoles, use of Eq. ([S]) has to be made when the "stationary" et 
and n are determined. This approximation is justified as we simply use the leading-order SPA and is valid in a regime 
where the orbit evolves slowly towards coalescence. It is, thus, an approximation around the point {no, eto) with the 
aforementioned requirements. 

These are the ideas so far. Here are the explicit equations up to IpN in conservative and only to leading-order in 
radiative dynamics, where n is unsealed (that means it has unit second^^). 



17 



n =M{n,et) = e° ^ — ^TtT^ ' (l^^) 



5(l-e?)V 

, 5 n7?et(l21e^ + 304) (Gm,n)5/3 
=^Ket) =-e ^5(^.^2) 5/2 • (193) 

In terms of combined quantities, related formulas for the simple Peters and Matthews approach have been published 
by Pierro and Pinto [13 and Appell's 2-variable hypergeometric function AppellFj^ has come to use. We will keep 
our own expressions for n and Cf and will express unsealed elapsed times as functions of the latter. Note that by an 
appropriate scaling the factor Gnit is re-absorbed in the time unit. 
The solution to dn/det ~ (dn/di)/(def/dt) reads 

/ _ 1) \ / N / i21e,o2 + 304^ 



with 640 and no as the value of et and n at the initial instant of time to, respectively. The elapsed time as a function 
of et and eto reads 



9519II8I/2299 



222173/2299co8/3^(Gmt)5/3 

4«/iq. „ /24 3 1181 43 o 121eto^ 
eto^^/i^AppellFi -; -, --— ; -; eto', 



19 2' 2299 19 304 

48/19, „ /24 3 1181 43 , 121ej2\] 
- ^^"^"AppellF, (l^; 2 ' j I ■ (195) 

18/19 /101 2 , Qf1^^ 1305/2299 

c„ (I21e,o + 304 no _ ^^^^^ 

(l-e?o) ^ 



A check will show that the right hand side has the dimension of time. 



C. Solution to the SPA condition equation, j 7^ 

It showed up that it is easier to solve Eq. (|186p for n instead of et and then to express e* — eto in terms of n — no. 
Expressed fully in terms of n, it reads 

= 27r/ - jn* ~ mk{n*,etin*))n* . (197) 
It is solved in two steps: first, solve the Newtonian and second: solve the IpN equation with the help of step 1, 

2tt f 

Step I : "n = -J- > (198) 

Step II : nJpN = 2nf - mHn*^,etin*^))n*^ ^^^^^ 

Step II drops out in case m = 0. The Newtonian "stationary eccentricity" for the frequency /, i.e. et(n^), can be 
found numerically or with the help of a perturbative solution scheme. It is a rather numerical issue to apply fixpoint- 
method-like iterative algorithms a la Danby and Burkhards' method [l^ to solve the Kepler equation, and a detailed 
analysis could indicate how many steps are necessary and reasonable towards the solution. Anyway, we promised to 
give analytical results and this will be done below. We are aware that there may exist better algorithms and refer 
the reader to the common literature of approximative solving methods. The function to be inverted for et — eto is the 
following, 



18 



= no+5^iHet-eto) + ?.9^'Het-eto)' + e~'g^'net-eto)^ + ... , (201) 
1 

9^'''' -T -5jr"("o,et,eto)|eto > (202) 

having introduced some smallness parameter e which will be set 1 after the calculation. The solution algorithm reads 
(defining k as the difference of and eto a-nd leaving out the "star") 



K 



(et-etor^, (203) 
«W = (204) 



5(1) 



-iy{(n-no)-?ff(2)(^[i])2| ^ (205) 

^^t^l = ^|(n-no)-f:e^-^(^)(^[^+i-''l)4 , (206) 

5 I p=2 J 

with some current solution order [N]. For convenience, we will give the first four orders of et — eto in terms of n — no- 
With the definitions 

/i := -1216^0 -183e2o + 304, (207) 

/2 376^0 + 29262^ + 96, (208) 



we obtain 



,[4] /„ „ \ /leto 



k'- ' ~ — (n — no) 



3/2 no 

2 /leto (370e^o + 344016^0 - 1318446^0 - 267206^0 - 56832) 



-e^ (n - no)^ j/iCto 

410700eti^ + 76370220ej^(^ 

+ (3257592723 - 48470/2) el^ ~ 3 (4124011/2 + 9090903688) e^^^ 
+168 (658637/2 + 276828486) e^^o + (9406764288 - 76265528/2) e^ 
+768 (131251/2 + 61327636) e^o + 18432 (267/2 + 494320) e^o 

-7274496 (7/2 - 1332)^ | 
x(l62/|n3)-^ 

-e^ (n - no)* < /leto 



+7597950006^4 



to 

+2119273605006^2 
-1850 (96940/2 - 10211818689) 6^° 
+ (459468074902815 - 6245O6868OO/2) 6^^ 
+90 (135716/1 - 42030994737/2 - 75142160154162) 6^^^ 
+8 (755339319/1 + 6763450991635/2 + 3192333977427390) e^o 
-3 (31474158721/1 + 56138837764240/2 + 8745121825435200) 6 
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+24 (7251390883/1 + 4696313592840/2 + 143898051631680) 6^^ 
-8 (35529727041/1 + 21737777348800/2 + 5492015496046080) 6^0 
+384 (94956913/1 - 19783158080/2 - 34649027828480) 6^ 
+30720 (2707169/2^ + 277826496/2 - 683225807616) e^^ 



.6 

-to 
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+ 18186240 (2319/1 + 594368/2 - 213546240) e^^ 
-1745879040 (35/| - 16576/2 + 1577088)^ | 

X {iMAf^Ji^y^ + 0{e^ {n - n^f) . (209) 

It is up to the reader to truncate this to some required order in em or to extend it in orders of e. Having found 
the "stationary" n\^-^ from Eq. (|199p and from e.1{n\^-^), one can obtain the associated t* by inserting this into Eq. 
(|195p . This in turn can be inserted into Eq. (|18ip to get the value of the phase at the stationary point t* . The reader 
should also note that the solution to Eq. (|186|) will introduce new IpN correction terms to the multipole moments of 



Eqs. ([TTaD - ([TTc]) . 

D. Solution to the SPA condition equation, j = and m 7^ 0. The pure periastron phase shift 

The stationary phase condition for the pure periastron-dependent terms (those with j — 0) reads 

$ = 27r/ — mkn 



27r/ - m7i(e*) 



0, (210) 



^g{et) 2^/-mn(e:)3 ^ \' ^0, (211) 

with n{et) taken from Eq. ()194p . Here, we proceed presenting all quantities expressed in terms of et- To find the 
solution to this equation analytically, we find it convenient to consider the perturbation algorithm from Danby & 
Burkhardt [20j to the fourth order. We need to have a nice initial guess for et, which we take from the first-order 
expansion of Eq. (|210p in et — e^o, 

[0] _ (242e[o + 1246^0 - 974e^o + 608eto) + 178e^o - 6696^0 - 784eto 

~ 3(19e4o- 284e?o-160) ' ^^^^^ 

noting that the case m = is excluded. This can be inserted into an iterative solution algorithm, which solves for 5 
in the expression 



(213) 

5(el°l+5) = 0. (214) 



This 5 is found with the help of the following procedure. 



<5i = -4, (215) 
9 

'52 = —tA^^ (216) 



\5i g 

'^3 = , , Ic , 1x2 III ' (^^'''^ 

5i = 1 , (218) 

9' + h5^9" + l5l9"' + ^A9"" 

We have e|^' = ej"^ + as the fourth-order solution to Eq. (|210p with quintic convergence, and again extract n(e['*'), 
t — to and so on. The case m ~ will be discussed below. 

E. The case j = and m = 0. Fourier transformation of a slow-in-time signal 

In Eq. (|177p . there is no fast oscillating term but only a slow variable of time to be Fourier transformed, 

r+oo 

Hf) = / A{t)e'''~^''^'Mt. (220) 
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The term A{t) depends on time only due to RR. Those terms are nontrivially dependent on time and have to be 
treated individuahy when they are requested analytically. In principle, one would have to express \E\, n and as 
explicit functions of time. That would include inversion of Appell functions or perturbation theory. 

However, for the case of inspiralling compact binaries, they will not be able to significantly contribute to frequencies 
in comparison to those with fast oscillating exponents as we compare typical time scales for one orbit and for the 
inspiral. We will therefore impose the following relation. 

'l(/)static ^ ^(/)stationary , (221) 

where /istatic means all Fourier integrals over terms where the mean anomaly - or equivalently, the time - does not 
appear explicitly. We state that the (j = 0,to = 0) Fourier domain terms almost vanish: 



/i(/)[,=o,™=o] « 0, (222) 

for the frequency domain of interest for the regarded detector. LISA for example, will hardly see those terms operating 
near no ~ O.OOlHz. Let us give an exemplary number to support this statement. The rate of change of the GW 
frequency / over one year will be Q 

1.6 X 10-9 f mt ( ri \ { U \^ f, , 73 ^ 37 4 



where fr is the radial frequency, given by fr = n(27r)-^. Let n ~ 10^'^, mi = 7712 = 1.4Mq and et = 0.1. Then, 
A/rr ~ 2xlO-^^IIz and the scaled energy loss is AEj^r ^ 3 x IQ-^^. 



F. The limit et 0: the quasi-circular case 

The quasi-circular inspiral has been discussed extensively in the literature, especially in Q which we have oftenly 
cited. For further information, see e.g. (2ll.[2^. In the limit et — > 0, all elements of our calculation simplify drastically. 
The following equations, 

M ^ M, (224) 
A{u) -> 1, (225) 
fiM ^ ^lUE), (226) 

gimiu) ^ girn{E), (227) 

(/.-0O ^ (l+k)M, (228) 

^-imrt> _^ g-im(/)Og-im(l+k)M (229) 

show that we have a simple prototype for the SPA for all the multipole moments. Eq. (|224p is Kepler's equation for 
quasi-circular orbits. The infinite summation series in Eq. (jl87p shrink to one single term, where the phase term has 
to be replaced by the one in Eq. (|229p . This is because the sin u terms will always have a factor and vanish in the 
case in question. The value of the phase and the angular velocity, the elapsed time and the resulting SPA integral 
can be taken one-to-one from lol. 



G. Some concluding remarks for the eccentric inspiral templates 

It is interesting to note how many parameters are included in the wave form. For the relatively simple case of 
circular inspiral, the templates used to have 

• tc, the time to coalescence, 

• (f>c, some phase instant, 

• Mc, the "chirp mass", 

• the symmetric mass ratio (rather important for higher- order corrections to the RR effects), and 

• i, the inclination angle of the orbital plane. 
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Because for eccentric orbits, both et and n will dictate the contribution to infinitely many frequencies already on the 
purely conservative level of EOM, both have to be regarded as parameters for the template. Thus, we have 

• TO, the total mass of the system, 

• 77, as before, 

• i as well, 

• no as the value of n at 

• eto as the value of et at Iq. 

Both no and eto will contribute to the time to coalescence, see (l^ for the value of what is called "lifetime". The 
parameter space has grown by one dimension, but the good news is that, for data analysis considerations, the ambiguity 
function is still maximized in view of (/)o in a considerably simple way (see (23j how to do this). 



VII. CONCLUSIONS AND OUTLOOK 



In this article, we provided the far-zone GW form, including IpN corrections to the orbital dynamics as well as to 
the amplitude. This was done by applying the IpN accurate QKP to the conservative dynamics first and to decompose 
each term in Fourier modes and in a second step to solve the Fourier integrals, modified by the leading-order effect of 
radiation reaction, with the help of the SPA method. The GW field is given in terms of tensor spherical harmonics 
and all terms are given in a purely analytical form, at least as they are solved to some required order. The inclusion 
of the fully analytically solved Kepler equation in terms of Bessel functions implies that single and double infinite 
summations appear. It is up to the user to restrict those summations to finite ones, as far as the required accuracy 
demands a minimal jmax and eccentricity expansion, due to the detector and other requirements in question. 

It is interesting to note that this approximation scheme is easily applicable to the case of spinning compact binaries 
with aligned spins and orbital angular momentum, including the leading-order spin-orbit interaction |24| . For non- 
aligned spins, the calculation of the Fourier domain will be structured more complicatedly, since the precession of the 
orbital plane will introduce another typical frequency in addition to the orbital, the periastron precession and the RR 
time scale frequencies. 

For a future publication it is intended to include higher-order RR terms to the "quadrupolar" contribution and, as 
it is highly demanded, to include the 2pN and 3pN point particle Hamiltonians into the dynamics. 
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Appendix A: Proof of some summation formulas 



In this appendix, we like to provide some proves of formulas we only listed in the previous sections. Throughout 
the remaining sections the eccentricity et we are using is simply called e. Let us start with the inverse scaled relative 
separation with some arbitrary positive integer exponent n. First, we perform a Taylor series expansion in e, 

1 ^ ^ (n-hTO- 1)! 

= 1+2^ TTT-^ r COS™ u 



(1 — ecosw)" ^ (n — l)\ to! 



, ^ (n + m-iy. e"' I 1 ^ fm\ „„A 



(n - 1)! to! \ 2™ ^ V ^ 

m=l ^ ' \ i=0 ^ 



and list the "factorial" function of the integer number n as 

(n -|- to — 1)! _ 



W{n + k-\). (A2) 
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To optically simplify this equation, we summarize the terms before cos ^ u with as follows, 

(„) (n + m — 1)1 1 e™ /to\ 

Pm,k ■= 1)1 '^\2^\k)' 

and write the sum with this definition: 



(A3) 



(1 — e cosu)' 



(A4) 



m=l fc=0 



To further markably reduce the complexity of this double sum, it is the task to find out which pairs of (m, k) lead to 
the same frequency ju and which have to be added to this frequency contribution: 



|m - 2fc| = j , 
=^ mi ~ 2k + j , 

A small table for the cos function argument may help (note that k < m always holds): 



(A5) 
(A6) 
(A7) 



m : k 





1 


2 


3 


4 


5 


1 


1 


-1 










2 


2 





-2 








3 


3 


1 


-1 


-3 






4 


4 


2 





-2 


-4 




5 


5 


3 


1 


-1 


-3 


-5 



The zero mode j = always appears at even numbers m, m = 2i with i = (1, 2, ...), and m — 2fc = 2i — 2fc = is 
satisfied by fc = i. Thus, we choose 



(A8) 



The other frequencies j > 0, appearing at m — j + 2i with i — (0, 1, 2, ...) lead to fc = i as well as fc = i + j, such that 



(A9) 



1=0 



Summarizing the elements for the first term, we have 

1 



(1 — ecosu)*^ 



6^"' cos JM 



(AlO) 



with the 6^"'' defined in Eqs. (|A8p and (|A9[) . Having this at hands, it is easy to compute the associated decomposition 
of 



sm u 



(1 — e cosw)'^ 



= I f + 6q"^ + ^ 5^"' cosju I smu 



^ OO 

(1 + 6^")) sinu + - ^ 5;."' (sin[(l - j>] + sin[(l + j>]) 



(All) 



This decomposition demands collecting the terms having the same frequency as well. Having m times u in the sin 
argument, following j will lead to sin(?T!,u) and — sin(mu) for sin(l — j)u in the first two lines and for sin(l + j)u in 
the third: 



m = (1 — j) — > J = 1 — TO , 

-TO = (1 - j) ^ j = 1 + TO, 



(A12) 
(A13) 
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j 


1 


2 


3 


4 


1+j 
1-j 


2 



3 
-1 


4 
-2 


5 
-3 



m = (1 + j) ^ j = m - 1 . (A14) 

Again, a small table may help: 

We clearly see that the first harmonic in this sum is only realized by j = 2 for the sin(l — j)u term, whereas j = 
will not contribute. Thus, 

c« oo oo 

b'j^^ (sin[(l — j)u\ + sin[(l + j)u\) — ^m^i sinmu — ^m+i sinmu — b'^^ sinu . (A15) 

j — l ni—2 ni—2 



Summarizing these terms, we write the Ihs of Eq. (|Alip as a simple sum. 



smu 



(1 — ecosu)" 



OO 


sin jM , 




2 2 


1 + 


- 

2 r^+i 


-!,<">) 



(A16) 



(A17) 



S^, := - (6^7, - b^a,) . (A18) 

Appendix B: Fourier representation of the sin mu and cos mu terms 

In the previous appendix, we have proven some formulas for terms appearing in the multipole expansion of the 
far-zone gravitational field. Now, we have to express simple series in sinmu and cosmu as trigonometric functions in 
ml. Using Eq. (fT5|) . we have 

^ oo oo / oo ' 

^ ^ rn—1 m—1 — 

c« / oo \ 

oc 

= 4")+^yl5")cosj7, (Bl) 

i=i 

4") := 1 + 6^"', (B2) 

-^S ■= fE^re') • (B3) 

\m=l / 

Similarly, getting help from Eq. (|16p . we obtain 

oo oo I oo \oo/oo \ 

a-TluY - E^^^"'--"--E^^"ME'^7'--^-0-E E^r^i^'M-^^^^ 

oo 

= Y^^f^^^jl. (B4) 
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Appendix C: Fourier representation of products of two sin series 
In section IIVI we provided a simple series representation of the term 



oo oo 



sin/sM ) ( Gra sinmM ) = X! X! ^4"'^™ (cos[fc - m]M - cos[fc + ?7i]M) 



\k=l 



k=l ■m=l 
oo 

jSS,[n] 



^P/^^1"1 cosjM. 

3=0 



(CI) 



Here we want to explain how the coefHcients p^^'^'I"! come up. First, we take a look at cos[fc — to]M. The zero mode 
appears only at fc = m. Thus, the absolute static part will be given by 



pSS..[n] C(n)r^ 



(C2) 



The higher harmonics j > will be realized hy k — m ~ j and k — m = —j and these equations solved for m give 
m = k — j (which only nonzero for k > j) and m = k + j . The first part, cos[fc — m]M, will give 



+'>fc <J-fe+j 



[l]pSS,[n] „ 1 ^ 

k=i \ 

The second part, — cos([/s + to]M), gives — cos jM for j ~ k + m, viz. m = j — k, such that it is equal to 



(C3) 



[2] pSSJn] _ { o(n)^ 



fc=i 



(C4) 



together giving 



pSS,[n] _ [l]pSS,l 
^j>0 - ^j>0 



^ oc 

[2] pSS,[n] _lsr 



k=l 



c(n)r' 1 c(")r' 



V 



(C5) 



k>3 



fc=l 



Additionally, we will provide the tables for part 1 and 2. 



k:m 


12 3 4 


k:m 


12 3 4 


1 


-1 -2 -3 


1 


2 3 4 5 


2 


10-1-2 


2 


3 4 5 6 


3 


2 10-1 


3 


4 5 6 7 


4 


3 2 10 


4 


5 6 18 


Numbers for part 1 


Numbers for part2 
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Appendix D: Fourier representation of products of a sin and a cos series 

What we did in the previous appendix is also necessary for the following term, 

Coo \/oo \_,oooo 

^^^"^cos/cMj I ^GfeSin/fcMj = - ^ ^ (sin[/c + m]M - sin[fc - m]M) 
/c=l / \*;=1 / fe=lm=l 

oo 

Pj ^ ' smjM . 



(Dl) 



Part 1, sin[fc + m\l will contribute to frequencies j = k + m for fc, to > 0. We already know the result from part 2 of 
the last section. 



(D2) 



k=l 



Part 2, — sin[fc — m\l will also have positive and negative contributions, for m = k — j (implying k > j) m ~ k + j, 
giving 



OO 

[2] pCS,[n] _ _i_ 
^3>l ~ 9 



fc=l 



\ for /c>j 



(D3) 



Summation of part 1 and 2 yields 



k=l 




for /c>j 



for j>l 



(D4) 



Appendix E: Appell's integral formula in the solution for the elapsed time 

The integral in Eq. p95p in section IVII can be solved with the help of the following integral representation of the 
AppellFl function (see http://dlmf.nist.gov/ or [25| for further information). 

^ \l-u[yHl-uvY^- ^ r(7) ^ AppellF, (a; A,/32;7;a:,y) . (El) 



Appendix F: Accuracy of finite sums 

In Sections IIIII and IIVI we provided decompositions of functions of some "elementary" type in terms of u which 
contain infinite summations. Naturally, for numerics it is important to know how many terms are needed to reach 
some desired accuracy. Considering compact binaries with small eccentricities only, one is allowed to expand the 
elementary expressions in powers of e and then to look how many terms are needed for the error to be shifted to 
0(e-^+i) with some finite M. 



1. Accuracy of basic elements 

We start with the basic definitions. The upper limits of the summation has to give a term of order ©(e^), thus 
the individual limit has to be matched appropriately, 

fiS.-^^^i:^~{":'\-0(e"-), (Fl) 



(n-1)! to!2'"V^' 
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M/2 



i=l 



i=0 



(7W-i)/2 



M 



Aiu)- 



(F2) 

(F3) 
(F4) 
(F5) 



In the last line we have used that the summation in Eq. (jF3p starts with i = and leaves no term if j > A4. The 
same quantity with sin it will also be truncated in the u domain, keeping in mind the definitions ([35]) and p6p and 
their dependency on the summation index j, 



smw 



gin) ^ 



1 + e 



r-(n) 

'i>i 



sinw 



/ M/2 \ (7W-2)/2 

y i=l J i=0 

^ '(7W-i+l)/2 (M-j-l)/2 

9 E (/5j-l+2i,i + /5j-l+2ij-l+i) ^ E \l^3 + l+2i,i + 1^, 



(F6) 
(F7) 



'i+i+2ij+i+j 



i=0 



M + l 



(F8) 
(F9) 

(FIO) 



For the Fourier representation, we remember Eq. I\37\i and take the expansion of the Bessel coefficients [I2[, 

^22fc+»fc!(fc + n)! ' 

for the determination of the hmit for our finite sums, 

777 

iT = - (J,-™(je)-J,+™0e)) = O(e^— ), 
J 

M+m 

^cosmw = ^ 7j"cosjM-0(e-^+i), 



M 



= (I + fe[,"')+^ylj"^cosjM + 0(e' 



(Fll) 

(FI2) 
(F13) 

(F14) 
(F15) 



In the third line above we used Eq. (jF5[) to truncate the number j in cosju. For the other relevant term with sinu 
we consider 



- -(J,_™(je)+J,+™(je))=0(e^-"), 

M+m 

sinmw = 5]"^ sin jM + 0(6-^+1), 



(FI6) 
(F17) 
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M + l 



ni—l 
M + l 



A{uy 



= 5]"^siiijM + 0(e-^+i). 



(F18) 
(F19) 



i=i 



In the third hne again, we used (jF9p for the index j in sin ju and in the last Eq. (jFlOp . We will also consider the 
expansion coefficients of v from Eq. (|23p . With the regular solution to Eq. ((24| at e = 0, 



1- 



their order is calculated to be 



Gm{e) = — Jm(me) +y^a^ [Jm_s(me) - Jm+s(me)] 

7D ^ ^ 



(F20) 



0(6™) . 



(F21) 



2. Accuracy of finite sum products 

The product of two sin series can be truncated taking into account the individual series' terms. In the following we 
will use the fact that those coefficients could already be truncated to some finite order Ai in the previous subsection, 
see Eqs. (|F14|) and (jF18|) . and that they are at least of some order of e themselves. 



4"^ sin fcM I ( ^ sin mM I = ^ p/^'I"' cos jM . 



(F22) 



J=0 



For clarity, we again write down the definitions and make use of the elaborated orders, where some abuse of notation 
is made, 







in) 
k 

(n) 



fc=l 



G, = 



(F23) 

(F24) 
(F25) 
(F26) 



fc=i 



k=l 



k=l 



jSS,l7l] 



fe=i 



for k>j 



k=l 



for j>l 



[fc=l k=j+l fc=l 

[fe=l k=j+l 
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M-j + l M+j + 1 

k=l k=j + l 

O(eJ-i). 



fe=i 



(F28) 



The individual final indices are evaluated when the maximal exponent of e reaches Ai. The last line is evaluated 
when one takes the smallest index k. The other product of cos and sin series is the following, 



A'k' cos /sM ^ Grn sm mM = ^ pj^'^-N sin jM . 



jCS,[n] ^ 

\k=l J \m=l / j=l 

The definitions are worked through immediately, again with some minor abuse of notation. 



(F29) 



P. 



CS\n\ 1 



J-1 



k+j 



k=l 



K for j>l 



k=l 



3-1 



k=j+i 



k=l 



3-1 



k=l 



for j>l 



fc=i+i 



k=l 
, for j>l 

i-1 



fc=l 



k=l 



k=j + l 



fc=l 



for J>1 



0(6^). 



Because of Eqs. (jF28P and ()F30|) . one can truncate the summations over the index j for 

1 . each ^j""* term at = Ai , 

2. each ^j""* term at j^vi = + 1, 

3. each p^^'^"^ term at Jm = M. + I, and 

4. each pj^"^''"' term at = A4. 
in Eqs. (|T55)) to PTO)) . 



(F30) 
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